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Abstract 

We discuss a new numerical method for the determination of ex- 



p cited states of a quantum system using a generalization of the 

Feynman-Kac formula. The method relies on introducing an 
ensemble of non-interacting identical systems with a fermionic 
statistics imposed on the systems as a whole, and on determining 
^ I the ground state of this fermionic ensemble by taking the large 

I time limit of the Euclidean kernel. Due to the exclusion principle, 

the ground state of an n-system ensemble is realized by the set of 
individual systems occupying successively the n lowest states, all 
of which can therefore be sampled in this way. To demonstrate 
how the method works, we consider a one-dimensional oscillator 
and a chain of harmonically coupled particles. 

Keywords: Feynman-Kac formula, antisymmetrization, excited states 
PACS: 02.70. Ss, 02.70.Tt, 05.10.Ln; 

Introduction 

A central object in the lattice formulation of field theory is the Euclidean 
kernel, which is obtained from the quantum amplitude by a Wick rotation. 
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For a (i-dimensional quantum field theory, the Euchdean kernel can be inter- 
preted as a partition function for a. {d + 1) -dimensional classical statistical 
theory. This theory can then be put on a lattice and simulated by Monte 
Carlo methods. 

The Euclidean kernel contains all the information about the quantum 
system, including in particular the excited states which in field theory cor- 
respond to particle states. Standard methods of extracting this information 
from numerical lattice data, commonly used for instance in QCD lattice 
spectroscopy, are based on the analysis of the large time behaviour of the 
Euclidean kernel |]T|. 

In this paper, we present a new method of determining excited states by 
studying the large time behaviour of the Euclidean kernel for an ensemble of 
identical systems with a fermionic statistics. The trick of using a fermionic 
ensemble of systems allows us to subtract the contribution from the ground 
state and focus directly on the contribution from the lowest excited state. 
The method can be recursively extended to the next excited states. 

The first part of the paper contains standard material, which we recall 
here mainly to introduce the notation for the second part of the paper. In 
this second part, we discuss the new method and present results obtained 
from applying it to the cases of a single harmonic oscillator and of a one- 
dimensional chain of coupled particles. 

We write most of our formulas in a quantum mechanical framework, which 
corresponds to a zero-dimensional {d = 0) quantum field theory. The gener- 
alization to [d > 0)-dimensional cases is straightforward and does not require 
a modification of the method or the underlying idea, although as we will see 
some practical limitations occur in this case. 

Standard methods 

The Euclidean kernel is defined as 



where (pniq) = are eigenvectors of the Hamiltonian H\n) = En\n) 

ordered hj En, Eq < Ei < . . . . The parameter r is real and non- negative. 
We set ^ = 1. The kernel can be expressed as a sum over paths propagating 
in time from an initial position at some initial time Tj to a final position 
qf at some later time Tf , t = Tf — Ti > 0, |^. Using the composition rule 





n=0 
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and applying it many times to equal time intervals e (A^e = r), one finds 



„ Af-l N 

Kr{qMi)= j n^^J W^Mk\<lk-l), (3) 

j=l k=l 

where Qf = Qn and = go- As it stands, the formula is of not much 
practical use, since the same unknown function K appears on both sides of 
the equation. The idea is now to use an approximation for small e, replacing 
on the right hand side K^^ with the semi-classical propagator /C^. As an 
example, consider a quantum mechanical system with the Hamiltonian H = 
+ V{q), which describes a particle in D dimensions. Then the semi- 
classical propagator /C^ is for small e given by 

}Ce{qk\qk-i) = (27re)-^/2 e-^^'^ (4) 

where 

AS, = U'-^^^^y + V{q,). (5) 



2 V e 

Using this, we can approximate the kernel (|^) for finite r = A^e by 

7V-1 N 



j=l k=l 



N-1 N-1 N 



K{qfh)^ I \{dqj \{lCMk\qk-i) = I Dq (6) 

where 

Dq = (27re)-^^/2 J] dq^ =C\{dq,, S = tY,ASk. (7) 

j=i j=i k=i 

For given A^ and e, the factor C is just a global normalization. In the limit 
N ^ oo with r = A^e = const, the formula (H) approaches the exact one (|]). 

The expression can be viewed as a partition function for a classical 
field {qj} on a one-dimensional lattice, j = 0, . . . , A^, weighted by the expo- 
nent of the Euclidean action : Dqe~^. One can simulate such an ensemble 
of {qj} using standard Monte Carlo techniques. It is convenient to impose 
periodic boundary conditions, i.e. to set qo = qN = q, and then to sum over 
q. The partition function for the ensemble of such periodic chains is : 

„ JV JV „ 

Z = Y[dqjY[lCMqk.i) = Dqe-'. (8) 

j=l k=l Aperiodic q 
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This function serves as an approximation of the quantum mechanical expres- 
sion 

I dq KMq) = jdq e"^^" { |0o(g) P + |0i(g)|V-(^-^«) + ...}. (9) 

On the one hand, the integrand dq K-r{q\q) corresponds to the probabihty 
distribution of positions q in the quantum system, but on the other hand 
it also describes the probability distribution of particle positions in the 
ensemble of classical chains. For each particle of the chain the distribution 
is the same because the partition function is translationally invariant. In 
the large r limit, the integrand dq — dq e~'^^''|0o('?)P is dominated 

by the contribution from the ground state. This is commonly known as 
the Feynman-Kac formula. The distribution of particle positions can be 
measured in Monte Carlo simulations of the classical chains. This can serve as 
a practical method for determining the ground state |(^o(q')P of the quantum 
system. 

In the same way, one can also measure the energy of the ground state. 
For each time slice k define an estimator of the classical energy : = 
Tk + Vfc, where Vk = V{qk) and Tk denote the potential and kinetic parts, 
respectively. The matrix elements of the squared momentum p"^ are given 
by (gfc+i — — therefore, the kinetic energy Tk depends in 

principle on three time slices. However, we can also introduce a one-slice 
operator for Tk by making use of the virial theorem, which tells us that on 
average {Tk) = \{qky {ik)) ■ The operator on the right hand side, qkV'ilk), 
is a one-slice operator and so, with these definitions, is now Ek- 

Denote the average over paths weighted by by (...). The average 
(Ek) does not depend on k, so we can average over slices E = Ylik to 
get an estimator with a smaller statistical error. In the large r limit, 

(E) = E, + {E^ - Eo)e--(^-^°) + . . . (10) 

approaches the ground state energy Eq. Since the average on the left hand 
side can be measured in the ensemble of chains, we thus have a method for 
computing Eq. 

In the MC simulations one is restricted to chains of finite length N . In 
order to increase r = A^e for finite N , one has to increase e. This intro- 
duces systematic errors since the deviations between and /C^ grow with e. 
Moreover, a small deviation from may accumulate and become amplified 
in the convolution of such factors in the time interval r = A^e (P). This 
error can be reduced by including higher order corrections in e to /C^ by 
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replacing the action with a sort of improved action : 

ASk AST' = ASk + ^e'iV\qk)f + ■■■ . (11) 

which can be obtained by introducing higher order corrections from the 
Baker-Campbell-Hausdorff formula in the semi-classical approximation (^. 
One can now make an optimal choice of and e so as to balance the sys- 
tematic and statistical errors. As it stands, the method may be applied only 
to the ground state since the contributions from higher states in (P,p!OD are 
exponentially suppressed. 

To determine the first excited state from the path integral approach, one 
has to remove the leading contribution. In the standard method this is done 
by calculating connected correlation functions for two different time slices ka 
and kb that are separated by a time interval At = u — Ta = e{kb — ka). For 
any one-slice operator O one measures the correlation function {0{Ta)0{Tb)) 
in the ensemble of chains, which is related to the following expression in the 
quantum system : 

{0{Ta)0{Tb)) = ^^e-("-^")^'"(m|0|n)e-^"^"(n|0|m). (12) 

m,n 

In the limit of large r, i.e. r ^ At, only the contribution from the ground 
state \m) = |0) survives in the sum. One gets in this limit : 

oo 

{0{Ta)0{Tb)) = 5^ |(0|O|n)|V^-(^"-^°) . (13) 

n=0 

The first term in the sum, |(0|O|0)p, is independent of At and can be sub- 
tracted by calculating the connected correlator : 

{{0{Ta)0{Tb))) = {0{Ta)0{Tb)) - {0{T^)){0{Tb)) . (14) 

For large At this is dominated by 

{{0{Ta)OiTb))) = |(0|O|l)|V^-(^l-^°) + . . . . (15) 

Thus, the difference Ei — Eq can be numerically determined by measuring 
the exponential fall-off of the correlation function (|15|) for large At. The 
freedom one has in choosing the operator O should be used in practice to 
maximize the coefficent |(0|(9|l)p relative to the coefficents for higher states 
n > 2. 

It is clear from the discussion above that this method requires a large 
separation of time scales -C At <^ r. This is a strong practical restriction. 
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Moreover, computations of connected correlation functions are in general 
much more time-consuming than calculations of averages. An additional 
difficulty comes from the fact that the signal-to-noise ratio is usually very 
small in the region of large At, where one has to carry out measurements of 
the exponential fall-off coefficient, Ei — Eq. 

To summarize, this method is much more demanding in computer power 
than the procedure for determining the energy of the ground state that we 
described before. 



Antisymmetrization over systems 

We now propose another way of subtracting the contribution from the ground 
state. We first consider an ensemble consisting of two non-interacting iden- 
tical systems, with a fermionic statistics imposed on the systems as a whole. 
The Pauli principle then forbids the two individual systems to be in the same 
state simultaneously. Thus, the ground state of the ensemble corresponds to 
one of the sub-systems being in its ground state and the other one occupying 
the lowest excited state. 

Denote the two systems by Q and R. They are independent and indis- 
tinguishable. Impose the fermionic statistics on them by defining the propa- 
gator for the twin system, K^'^\ as the antisymmetrization of the individual 
propagators : 

Kl'\qj,rj\q,,n) = ^ {KAqj\q,)KArf\n) - KArfh)KAqf\u)} • (16) 

It follows from (|[) that this function fulfills the same kind of composition 
rule : 

Ki%{qf,rf\qi,ri) = j dqdrKl^\qf,rf\q,r)K^J\q,r\qi,ri). (17) 

Thus, we can repeat the same approximation scheme as for the individual 
propagators (H). We get 

N-l N 

Kf\qf,rf\qi,ri) ^ / JJ rfg^c/rj JJ £f ^(gfc, r^lgfe-i, rfc_i) , (18) 

A — 1 1 1 



3=1 k=l 



where 

lCf\qk,rk\qk-i,rk-i) = ^{lCe{qk\qk~i)JCe{rk\rk-i)-JCeirk\qk~i)JCe{qk\rk-i)} ■ 



(19) 
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Now we can define the partition function for tfie QR ensemble with periodic 
boundary conditions, analogously to The partition function is related to 
the underlying quantum mechanical quantities as follows : 

Z^^-* = J dqdr Ki^\q,r\q,r) 

= Jdqdr e--^o''||$(2)(^^^)|2 ^ |^(2)(^^ ^^^^-.(^f)-^'^') ^ _ _ _ | 

^ e-^^o^ Jdqdr \¥o\q,r)\^ 

(20) 

where now ^^^\q, r) are wave functions of the twin system and E^^ are the 
corresponding energy levels. The wave functions are constructed as Slater 
determinants of the wave functions of the individual systems. In particular, 
the ground state is $q {q,r) = {0o('?)0i(^) ~ 0i(Q')0o(^)}! and its energy 

is given by Eq = Eq + Ei. This provides us with a practical method of 
determining the lowest excited state |0i(q')P and the corresponding energy 
El. The probability distribution of positions in one of the two systems, 
P^'^\q), can be expressed in terms of the wave functions of the individual 
system : 

p(2)(g) = ldr |$?^(g,r)r = 1{\Mq)\' + \Mq)\'} (21) 
Both P^^\q) and P^^\q) = |0o(g)P 

can be measured numerically using the 
MC method, and we can put them together to find 

|0i(g)P = 2P(2)(g)_p(i)(g) . (22) 

Similarly, by subtracting the ground state energy Eq of an individual system 
from the ground state energy Eq of the twin system we can determine the 
energy Ei = Eq — Eq of the first excited state of the individual system. 

This method can be extended recursively to determine higher excited 
states as well. Namely, the energy Ek of the k-th excited state can be com- 
puted as Ek = Eq'^'^^^ — E^^ , the difference of the ground state energies of the 
ensembles composed of k and k — 1 antisymmetrized copies of the individual 
system, respectively. Similarly, for the /c-th excited state, 

= {k + l)P^'+'\q) - kP^'\q) . (23) 
where P^^-* is an appropriately normalized probability distribution 

pW(g)= [ \^^\q,r,s,...)\^ drds ... (24) 
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Figure 1: The histograms show the numerical results for P2{<l) and Pi{q) 
from the MC simulations. The solid line represents the theoretical curve 

for the ensemble of k antisymmetrized copies. In the Monte Carlo simulations 
of the ensemble of k copies, the probability distribution P^''\q) is measured 
straightforwardly as a probability distribution of particle positions q in one 
of the k copies. 

Results 

Let us first illustrate how this method works for the one-dimensional har- 
monic oscillator : V{q) = q'^/2. In figure we show the probability distri- 
butions P^^\q) = |0o(q')P and P^'^\q) in simulations with a single system 
and with the ensemble of two antisymmetrized copies, for N = 128 and 
e = 0.0625 (corresponding to r = 8). The solid line going through the data 
points of P^^\q) is given by the ground state function of the oscillator : 
I'/'ol'?)!^ = ^^^ '^ • '^^^ difference of the two numerically obtained curves 
2p(2)(g) _ p(i)(g) is shown in fi gure 0, where it is compared to the function 
= ' which describes the first excited state squared of the 

oscillator. The agreement is perfect. 

Physically, the appearance of the valley in the P^'^\q) distribution re- 
sults from the fermionic repulsion between the two systems. The difference 
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Figure 2: The histogram shows the normahzed difference 2P2{q) — Pi{q) 
of the data from figure |1|. The solid line represents the theoretical curve 

\4'i{q)\'^ = 2P^^^(g) — P^^\q) leads to a zero value of the probability. 

The measured values of the energy of the ground state are Eq = 0.501(3) 
and = 2.006(6), which gives El = 1.505(7) in agreement with E^ = n+^. 

The method also works very well for higher excited states. As an example 
we show, in figure ^ the probability distributions P^^\q) and P^^\q) and, in 
figure^, the difference AP^'^^q) — 3P^^\q). The resulting curve is compared 
to the function |03(g)P = 2^! ^-^3 (^)'^~^ ' "^^ere H^lq) = 8q^ — 12q is the 
third Hermite polynomial, drawn as a solid line. Again, the agreement is 
very good. 

As a second example, consider now a one-dimensional chain of harmoni- 
cally coupled particles with the Hamiltonian 

(25) 

The index j runs over positions in the chain. The harmonic coupling constant 
between neighbours is l/e"^. As before, using the path integral formulation 
we first write the kernel as a convolution of N propagators (^, and then 
approximate JC^ in the limit e — *^ 0. Similarly to (P), we get a classical 

partition function for the ensemble of fields qkj, distributed now on a two- 
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Figure 3: Numerical results for ^4(5) and Psiq)- 
dimensional lattice and weighted by the action 

Se = J: \ [ '"^'■\- ""■' )' + \ + V(,,.,) (26) 

with the integration measure ~ Yl dQk,j (0); where the index k runs over time 
slices. Note that despite the apparent similarity, the two parameters e and 
e have a different meaning : the former is a time-slicing parameter, whereas 
the latter is a harmonic constant of the chain. 

One way of looking at the chain is from a quantum mechanical perspec- 
tive. In this case one fixes the harmonic constant e and for given r makes 
the time-slicing dense, N ^ 00 and e = r/iV — >■ 0. Alternatively, one can 
adopt a field-theoretical viewpoint by re-interpreting the parameter e as a 
spatial lattice constant and keeping it proportional to e : c = e/e = const. 
In this case, the pair {k,j) can be viewed as a space-time index. In the naive 
continuum limit e — > 0, the two terms in the sum (^) become a Lorentz (ro- 
tationally) invariant expression {dtqY + c^ld^q)"^, and therefore so does the 
whole theory. Thus, the expression ( ^61) can be regarded as a discretization 
of the Euclidean action of a field theory, and the partition function as the 
Euclidean kernel of a quantum field theory. We see that the prescription for 
the quantization of the theory is almost the same as in quantum mechanics. 
It is commonly treated as a non-perturbative definition of quantum field 
theory. 
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Figure 4: The normalized difference 4P4(g) — 3P3(g) of the data from figure 
^. The solid line represents the theoretical curve |03(Q')p. 

A difference between quantum mechanics and field theory appears, how- 
ever, when one wants to introduce the e corrections. The formula (|ll]), which 
works for quantum mechanics, does not work for quantum field theory since 
it breaks the symmetry between space and time. For quantum field theory 
one therefore has to use a different scheme (an 'improved action') to intro- 
duce lattice spacing corrections to the action which preserve the Lorentz 
invariance of the underlying continuum theory [^. 

Here we will consider only the simplest case, namely a quantum chain 
with the potential V{q) = which we simulate on a periodic (1 + 1)- 

dimensional lattice with nodes in the temporal direction and n in the 
spatial one. To each node we assign an appropriate particle displacement q. 

The spectrum of this simple Hamiltonian is known. For convenience, 
we will write formulas only for odd n = 2 A + 1. The Hamiltonian can 
be diagonalized in the Fourier modes Qa, a = 0,1, . . . ,A. The modes for 
a > are twice degenerate, each having a left mover and a right mover. The 
frequencies of the modes are 



= \/ 1 + ^ sin^ — , a = l,...A. (27) 
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n = 1 


2 


3 


5 


7 


= 16 


0.981(3) 




0.967(2) 






32 


1.001(3) 




0.993(4) 






64 


1.001(3) 


0.998(2) 


1.001(7) 


1.00(5) 


1.04(7) 


128 


0.995(3) 




1.02(3) 







Table 1: The mass gap m for various N and n {t = 5, e = 1). 



This leads to a ground state (vacuum) energy 

1 ^ 

Eo = - + ^^„ (28) 

and the following energies of one-particle states, numbered by the Fourier 
mode index a, 

AEi{a) = Ei{a) - Eo = , (29) 

where the particle mass is m = AEi{0) = 1. 

For a -C n, the equation (|29|) reproduces the standard forumula AEi = 
a/ + with the momentum p = 27ra/{en). 

The mass m can be calculated as the difference between the ground state 
energy of the fermionic twin system and that of one individual system. In 
practice, for finite r such an estimator also includes contributions from non- 
zero momentum states. Although they vanish exponentially with r, for finite 
r they may systematically lead to results that are somewhat overestimated.^ 
To minimize this unwanted effect one should concentrate the analysis only 
on the zero momentum mode Q = Qa=Q- The effective action density for the 
zero-momentum mode is given by (^, with q replaced by the zero mode Q. 

To determine the energy of Q, we simulate the whole (1 + l)-dimensional 
system, estimate for finite r the probability that the system is in the zero 
momentum mode, and measure quantities related to Q only. The results are 
summarized in tables 1-2. 

Table 1 contains the results for the particle mass. As we can see, the 
results agree with the theory; in other words, the method works. Unfortu- 
nately, the method is limited by the sign problem, which enters the game 
as a result of the antisymmetrization (p!9D, causing the integrand of ([T8|) to 

^The effect grows with n because the energy separation of momentum modes decreases 
as n increases. For example, at t = 5 the maximal effect we measured was of order 10% 
for n ~ 1 . 
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n = 1 


2 


3 


5 


7 


= 16 


1 




0.2221(3) 






32 


1 




0.0912(5) 






64 


1 


0.2044(6) 


0.0436(4) 


0.0120(3) 


0.0104(3) 


128 


1 




0.026(1) 







Table 2: The average sign for various and n ( r = 5, e = 1). 



not be positive in general. It forces us to first change the MC weights by 
taking their absolute value |/C| of the integrand in ([ISD , use these modified 
weights in the simulations, and then eventually to include the sign of /C in 
the estimators of the measured quantities : 

(C sgn(£))|^i 

(O) = ^ ^ 1 ^'1^1 . (30) 

(sgn (/C))|£| 

The sole exception to this is the quantum system with only one degree of 
freedom, for which it is possible to show that the sign of the weights is always 
positive due to cancellations appearing in the two-step transfer matrix. 

Generically, one expects the sign problem to create an exponential growth 
in the computer demands both in n and A^, because the computer time 
required is roughly inversely proportional to the average sign, which typically 
approaches zero exponentially. To test this, we measured the average sign as 
a function of n and A^. As one can see from table 2, it does indeed decrease 
with growing A^. On the other hand, as can be also seen in the table, for 
small A^ the sign is far from zero, but then the results for the particle mass are 
biased by systematic errors. It is a practical question whether there exists 
a window in A^ that has both small enough systematic errors and a large 
enough average sign. This might be the case for systems which have excited 
states with relatively large masses that decay over a small number of time 
slices, such as for instance SU (3) . 

Conclusions 

In conclusion, we presented a new method for the determination of excited 
states in quantum mechanical systems. For quantum mechanics with one 
degree of freedom, the method allows to determine recursively the lowest 
states and their energy. The method can be extended to systems with more 
degrees of freedom, but in doing so one encounters the sign problem which 
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in many cases can create practical limitations. It may be that applying the 
recent ideas for defeating the sign problem [0 to the method presented here 
would lead to an efficient method for field theoretical applications as well. 
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